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ABSTRACT 

We study the force of dynamical friction acting on a gravitating point mass that 
travels through an extended, isothermal gas. This force is well established in the hy- 
personic limit, but remains less understood in the subsonic regime. Using perturbation 
theory, we analyze the changes in gas velocity and density far from the mass. We show 
analytically that the steady-state friction force is M V, where M is the mass accretion 
rate onto an object moving at speed V. It follows that the speed of an object experienc- 
ing no other forces declines as the inverse square of its mass. Using a modified version 
of the classic Bondi-Hoyle interpolation formula for M as a function of V, we derive an 
analytic expression for the friction force. This expression also holds when mass accretion 
is thwarted, e.g. by a wind, as long as the wind-cloud interaction is sufficiently confined 
spatially. Our result should find application in a number of astrophysical settings, such 
as the motion of galaxies through intracluster gas. 

Subject headings: hydrodynamics — ISM: general — galaxies: kinematics and dynam- 
ics — stars: kinematics 

1. Introduction 

A gravitating mass that traverses a sea of other particles builds up an overdense wake behind 
it. This wake tugs back on the mass, providing an effective drag. The background sea may itself 
consist of non-interacting point masses. For this collisionless case, Chandrasekhar (1943) first 
derived the dynamical friction force. His celebrated result has since found application in a great 
many astrophysical problems, ranging from mass segregation in dense star clusters (Portegies- 
Zwart & McMillan 2002) to planet migration through interaction with planetesimals (Del Popolo, 
YesUyurt, & Ercan 2003). 

The background environment may also be an extended gas cloud. This type of dynamical 
friction has also been invoked in a variety of contexts, such as black hole mergers in galactic nuclei 
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(Dotti, Colpi, & Haardt 2006), the heating of the intracluster medium by infalling galaxies (El- 
Zant, Kim, & Kamionkowski 2004), and the migration of giant planets within a protoplanetary disk 
(Ogihara, Duncan, & Ida 2010). When the ambient medium is a gas, pressure gradients influence 
the formation of the wake behind the gravitating object. Surprisingly, the general determination 
of gaseous dynamical friction, for arbitrary Mach number of the perturbing mass, has not yet been 
achieved. There is substantial agreement when the mass is traveling hypersonically with respect to 
the gas. In this limit, the force varies as V~ 2 , where V is the speed of the perturber (Dokuchaev 
1964; Ruderman & Spiegel 1971; Rephaeli & Salpeter 1980; Ostriker 1999). 

In all studies thus far, the authors first calculated the properties of the wake by treating it 
as a linear perturbation of the background gas. The density perturbation is symmetric upstream 
and downstream when the mass is moving subsonically (Dokuchaev 1964), leading Rephaeli & 
Salpeter (1980) to conclude that the friction force is zero in this case. Ostriker (1999) obtained 
the force through direct integration over surrounding fluid elements, using their respective density 
enhancements. She added the constraint that the projectile's gravitational field only be switched on 
for a finite time interval At. With this device, she first found a nonzero result even in the subsonic 
regime. The force increases with V, and logarithmically diverges at a Mach number of unity. 

Interestingly, the quantity At does not appear in Ostriker's final expression for the subsonic 
force. This fact indicates that the artifice of a finite time interval was unnecessary and that a 
steady-state analysis is applicable. Indeed, the force attains a steady-state value in the numerical 
simulations of Sanchez-Salcedo & Brandenburg (1999). The divergence at a Mach number of 
unity in the analytical expression further suggests that physical understanding of the problem is 
incomplete. 

In this paper, we revisit the subject of dynamical friction, concentrating entirely on the less 
studied subsonic case. We take the perturbing body to be a point mass M traveling through an 
initially uniform gas. The previous studies cited also ostensibly dealt with point masses, in the 
sense that the physical size of the body was ignored. However, it was assumed, either tacitly or 
explicitly, that the object's radius R far exceeds the accretion radius r acc , conventionally defined as 
r acc = 2G M/V 2 . It is true that when R S> r acc , the gravitational force from the object is so weak 
that mass accretion by infall is negligible. Under these circumstances, however, the primary drag 
on the body is not from dynamical friction, but from direct impact by the gas, a fact sometimes 
overlooked. 1 

The conventionally assumed inequality marginally holds in one situation commonly envisioned, 
galaxies within intracluster gas (R ~ r acc ~ 10 24 cm). However, it fails badly in other contexts, 
e.g., supermassive black holes within galaxies (R ~ 10 10 19 cm) or gas giant planets 



1 Ruderman & Spiegel (1971) recognized that both drag forces act on galaxies moving supersonically through 
intracluster gas (see their eq. 5). They extended the linear analysis of the flow into the nonlinear regime, utilizing a 
similarity solution. However, their focus was the X-ray emission from the wake and bowshock, rather than the actual 
motion of the galaxies. 



-3 - 



inside circumstellar disks (R ~ 10 10 cm, r acc ~ 10 13 cm). When R <C r acc , as we assume here, 
dynamical friction is indeed the main drag force. The relative density enhancement in the wake is 
not small, as needed for linear theory (see, e.g., Kim &: Kim 2009), and mass accretion cannot be 
neglected. 

Our analysis indeed pivots on the fact that the transfer of linear momentum from the back- 
ground gas to the object, which underlies the friction force, is closely related to the transfer of 
mass. The problem of gas accretion onto a moving body was addressed in a classic series of papers 
by Hoyle & Lyttleton (1939), Bondi & Hoyle (1944), and Bondi (1952). The final result for the ac- 
cretion rate, applicable for all Mach numbers, is the interpolation formula offered by Bondi (1952). 
While not derived rigorously, the formula matches known results in the hypersonic and station- 
ary limits, and is broadly consistent with numerical simulations (see Ruffert 1996, and references 
therein). 

The strategy in our paper is to determine, using perturbation theory, the density and velocity 
of the gas. However, we focus not on the wake, as in previous studies, but on a region far from the 
object, where its gravity is relatively weak. Extending the perturbation analysis into the nonlinear 
regime, we calculate the net momentum flux onto the accreting object and derive analytically that 
the force from dynamical friction is MV, where M is the mass accretion rate onto the object. 
Adopting an analytic form for this rate, the drag force follows. This force first rises with V and 
then falls, remaining finite at all Mach numbers. Moreover, there is a contribution from the direct 
accretion of fluid momentum onto the body. This contribution is absent in the stellar dynamical 
problem, but is here comparable to the gravitational tug from the wake. 

In Section 2 below, we introduce a perturbative series expansion to analyze the far-field density 
and velocity. In Section 3, we use this expansion to derive a hierarchy of dynamical equations, of 
which we need only solve the first two sets. Section 4 shows how the mass accretion rate is connected 
to solutions of our second-order equations. In Section 5, we similarly relate the friction force to 
these solutions, and derive the central connection between this force and the accretion rate. Using a 
modified version of the Bondi interpolation formula for the latter, we find explicitly the deceleration 
of an isolated mass in Section 6. Finally, Section 7 compares our result with existing numerical 
simulations and suggests future avenues of inquiry. 

2. Outer Flow: Method of Solution 

2.1. Physical Assumptions 

Let the gravitating mass M travel in a straight line with speed V through the extended gas 
cloud. Following previous analytic studies of dynamical friction, we assume the gas to be isothermal, 
with an associated sound speed c s . Very far from the mass, the density is spatially uniform and has 
the value pq. We choose a reference frame whose origin is anchored on the perturbing mass. In this 
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frame, it is the gas that has speed V far from the mass. We let the gas velocity be directed along 
the z-axis, and employ spherical coordinates r and 9 (see Fig. 1). We now seek a steady-state, 
axisymmetric solution for the flow, which is taken to be inviscid. We neglect the self-gravity of the 
gas, and assume that each fluid element feels only a pressure gradient and the gravitational pull of 
the point mass. 

Strictly speaking, there is no steady flow, as this mass decelerates and V continually changes 
(e.g., Fathi 2010). What, then, is the meaning of the force we are calculating? Imagine the object 
being dragged by a massless string through the gas at the fixed speed V. After a long time, a 
steady-state flow is indeed established throughout the surrounding gas, and the tension in the 
string approaches a constant value. This limiting tension is the dynamical friction force being 
calculated here. 

Return now to the actual case, in which there is no string and the mass decelerates. As stated 
previously, there is no global, steady-state flow. The flow is quasi-steady, however, within some 
distance over which the altered motion of the mass is communicated by sound waves. We shall 
quantify this distance later, after we have established the flow assuming steady-state conditions. 

One important property of the flow is that it is irrotational. Euler's equation in steady state 
may be written 

u x u = VB , (1) 
where u is the fluid velocity, u = V x u is the vorticity, and the Bernoulli function B is 

B^-u 2 + c In . 2 

2 \poJ r 

Both the fluid speed and density approach constant values far from the mass. Hence, B is a spatial 
constant throughout the flow, and 

u x u = . (3) 

Since it is a poloidal vector, the vorticity u is toroidal. The last equation then implies that u> = 0, 
as claimed. We will not need to invoke the irrotational character of the flow until Section 5, when 
we explicitly evaluate the dynamical friction force. 

Throughout our analysis, it will be more convenient to employ, not the vector fluid velocity 
u(r,9), but the scalar stream function ip(r,0). We may recover the individual velocity components 
from the stream function through the standard relations 

1 d± (A s 

Ur ~ pr* sine ae ' K ) 

-1 dt/j 

ue = i—R ~o ' (5) 

p r sin t) or 



where p = p(r, 9) is the mass density. The velocity, as given by equations (4) and (5), automatically 
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Fig. 1. — Spherical coordinate system centered on a gravitating body of mass M. The gas is 
isothermal, and its velocity far upstream is /3c s (/3 < 1). The upstream direction corresponds to 
6 = 7r, and downstream to 9 = 0. While the figure shows the mass to have a finite physical size, 
we assume it to be a point particle in our analysis. 



- 6 - 



obeys mass continuity: 



= V-Gou 

i a , 2 , id 

— (pr u r ) H — - 

or r sin ft ov 



= Z2^z(p r2u r) + Z^Z^l^(P s[ndu e) ■ ( 6 ) 



2.2. Perturbation Expansion 



Far from the mass, as the density approaches po, the velocity has only a z-component, which is 
V. Equivalently, we have in this limit u r ~ V cosO and uq ~ — V sin#. It follows that the far-field 
limit of the stream function is 

po V r 2 sin 2 6 



(7) 



For a more complete analysis of the flow in this region, we take equation (7) to represent the 
leading term of a perturbation expansion. Introducing the sonic radius r s = GM/c 2 s , we first 
rewrite equation (7) as 



where /3 is the Mach number of the projectile mass: 



2 • 2 

r \ sin i 



(8) 



(9) 



Our perturbation expansion of ip is then given by 



Po c s r s 



where 



h - + h - + /o + /-i - 



/2 



+ 



/3 sin 2 



(10) 



(11) 



and where /i, fo, etc. are still unknown, nondimensional functions of /3 and 9. Similarly, we 
expand the density as 



Po 



1 + 9-1 - 



+ 5-2 — + 5-3 — 



+ 



(12) 



Here, <?_i, g_2, <?-3, etc. are also nondimensional functions of (3 and 0, all yet to be found. 
Both expansions are only valid for r 3> r s , the inequality that defines our outer region. We further 
assume that the physical radius of the object obeys R <S r s . Since the motion is subsonic (V < c s ), 
it follows that R r acc , so that mass accretion is significant. 

At this point, it is convenient to cast all our variables into nondimensional form. We let the 
fiducial radius, density, and speed be r s , po, and c s , respectively. Similarly, the stream function 
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is normalized to poc s r1. Then equations (4) and (5) relating the velocity to the stream function 
remain the same nondimensionally. We shall not employ a new notation for nondimensional vari- 
ables, but make it clear whenever we switch back to dimensional relations. With this convention, 
the nondimensional expansions for the stream function and density simplify to 



By adopting these perturbation expansions, we have effectively limited our analysis to the 
subsonic regime. For /3 > 1, we expect the fluid variables or their derivatives to be discontinuous 
across the Mach cone, whose opening angle is given by sm.9 = j3~ l (see, e.g., Ruderman & Spiegel 
1971). It would thus be necessary to adopt two separate expansions for ip and p, one inside and 
one outside the Mach cone. To avoid this complication, and since we are primarily interested in 
the subsonic regime in any event, we assume that ft < 1 and retain the single expansions. 



By symmetry, the upstream axis of the flow, defined by 9 = tt, is a streamline for any /?- value. 
That is, ip(r, tt) is independent of r. The actual value of ip(r, tt) is immaterial, reflecting the fact that 
the full function ip(r, 9) can have an arbitrary additive constant without affecting the velocities. For 
convenience, we set tp(r, tt) = 0, and note from equation (11) that /2(tt) already vanishes. From 
equation (13) for the general expansion, we require that fi(ir) = 0, for i = 1,0, —1, —2, etc. 

A second set of boundary conditions pertains to the behavior of the velocity u. Let us focus 
again on the upstream axis. The righthand sides of both equations (4) and (5) contain sin 9 in the 
denominator. Since the density p is finite at 9 = tt, where sin# vanishes, both dip/ 89 and dip/dr 
must tend to zero as 9 approaches tt, at least as fast as sin#. 

Considering first dip/dO, we see that df2/d0 = j3 sin 9 cos 9, which properly vanishes. We must 
further demand that f- (tt) = 0, for % = 1,0, —1, —2, etc. Turning to dip/dr, the term involving 
fi still goes to zero, while fo(9) itself vanishes when taking the r-derivative of ip. We are already 
requiring that fi (tt) = for all other i. Thus, the stipulation that ip (r, tt) = implies that dip/dr 
also vanishes at 9 = tt, so that ug does not diverge. 

Approaching the downstream axis, sin 9 again vanishes as 9 goes to zero. By analogous reason- 
ing, we require that f[ (0) = for i = 1,0, —1, —2, etc. To ensure the regularity of uq, we further 
need fi (0) = 0, for i = 1, —1, —2, etc. Again, the term disappears when taking the r-derivative 
of ip, and there is no a priori restriction on /o (0). Indeed, this quantity sets the mass accretion 
rate onto the moving body, as we later demonstrate. In summary, our boundary conditions are: 



fi (tt) = // (tt) = f! (0) = 0, for i = 1, , -1, -2, etc., and /; (0) = for i = 1, -1, -2, etc. 



^ = h r 2 + fi r + f + /_! r- 1 + ... , 
p = l + 5-l r' 1 + 9-1 r ~ 2 + 5-3 r~ 3 + 



(13) 
(14) 



2.3. Boundary Conditions 
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3. Outer Flow: Results 

3.1. First-Order Equations 

Our inviscid flow obeys Euler's equation, which we write in spherical coordinates. The r- and 
^-components of this vector equation are 



^ dUr_ + Ug^du^ _ V% _ _\df>_ \_ 

dr r 89 r p dr r2 



dug ug dug u r u e 1 dp 

Ur ~5 - + 55" + = 55 ' ( 16 ) 

dr r do r pr ov 

where u r and ug are given in terms of ip by equations (4) and (5). Our strategy is to substitute the 
perturbation expansions (13) and (14) into Euler's equations. For each power of r, we demand that 
its coefficients match. In this way, we will obtain a hierarchy of coupled equations for the functions 
/ and g. 

Before proceeding, we first note that u r and ug are both proportional to p~ l . To avoid expand- 
ing inverse powers of the density, we multiply equations (15) and (16) through by p 3 . Replacing the 
velocity components by derivatives of ip results in complex expressions that we shall not write out 
in full. We simply note, as an example, that the first lefthand term in the r-component of Euler's 
equation is 



du 



r 



p fd^\ 2 1 dp fdi/)\ 2 1 d^j d 2 ^ 



dr r 5 sm 2 \ d6 J r 4 sin 2 dr \ d8 J + r 4 sin 2 6> d9 drd8 ' ^ 

After substituting the series expansions for ip and p, we find that the highest power of r is r~ 1 . In 
this case, all the coefficients on both sides of Euler's equations vanish identically. 

Matching the coefficients of r -2 , we obtain the first-order equations. From the r-component 
of Euler's equation, we find 

-Pf" ~ Mi + /3 2 sin#cos0 g'_ x + (/3 2 cos 2 # - l) g_ x + 1 = , (18) 

while the (^-component yields 

(1 - (3 2 sm 2 9) g'_ x - /3 2 sin0cos# 5 _i = . (19) 

In both of these equations and those that follow, a prime denotes a ^-derivative. 

These equations govern the first non-trivial terms in the expansions for ip and p. Their solu- 
tion, therefore, must be equivalent to that obtained through the more traditional, linear analysis. 
Integrating equation (19), we find 

9 -' = h R2- 2 Ml/ 2 ' (20) 

(1 - /3 2 sur#) 1 
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where C is a constant, as yet undetermined. In the subsonic regime, the denominator on the 
righthand side does not vanish for any 9, and g-\ remains finite. 

Substituting this expression for g_i into (18) gives the equation obeyed by f±: 

f" + f - 1 ° (1 " (21) 

A particular solution of this equation may be found through the method of variation of parameters. 
Adding the two homogeneous solutions yields 

1 C (1 - /3 2 sin 2 fl) 1/2 , 
fl = p + £>cosfl + Esin9 , (22) 

where D and E are also constants. 

We proceed to evaluate the constants through application of the boundary conditions. The 
requirement that /i(vr) = gives 

C = 1 - PD . (23) 

Similarly, we have /i(0) = 0, yielding 

C = 1 + PD . (24) 

It follows, from these last two relations, that C = 1 and D = 0. Finally, we have /{(vr) = 0, from 
which we infer that E = 0. It may be verified that the boundary condition /{(0) = is then also 
satisfied. 

In summary, the first-order density and stream function perturbations are 

9 ~ X = h J' 2 /A 1/2 ' (25) 

(1 - (3 2 sm z 8) ' 

l _ (i _ /?W0) 1/2 

fl = ~ p — ■ (26) 

Notice that g'_ 1 = at both 6 = and tt, implying that the density profile is flat (i.e., does not 
have a cusp) on either the upstream or downstream axis. Our expression for g_i is consistent 
with the linear density perturbation obtained by Dokuchaev (1964), Ruderman & Spiegel (1971), 
and Ostriker (1999). Notice that g_i(n/2) diverges as /3 approaches unity, signifying the birth of 
the Mach cone. Our f±, in combination with gr_i, yields the linear velocity components given in 
equations (25) and (26) of Dokuchaev (1964). 

Figure 2 displays the streamlines (solid curves) and isodensity contours (dashed curves) for the 
outer flow, including only the first-order perturbations. The light, dotted circle marks the sonic 
radius, r = 1; the solution is only accurate well outside this sphere. Notice how all curves and 
contours are symmetric about the 9 = tt/2 plane. The streamlines, in particular, show the fluid 
veering toward the mass, but then turning away again. We cannot detect true accretion of mass or 
linear momentum until we include the next higher-order perturbations. 
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Fig. 2. — Streamlines (solid) and density contours (dashed) for the /3 = 0.5 flow, including only first- 
order perturbations. All quantities shown are nondimensional. The density contours correspond 
to p = 1.2, 1.4, and 1.6, while adjacent streamlines enclose equal mass fluxes. The inner dotted 
circle has the sonic radius. Since the streamlines and density contours are symmetric upstream 
and downstream, we cannot determine the true accretion rates of mass and momentum without 
higher-order perturbations. 
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3.2. Second-Order Equations 

We next equate coefficients of r~ 3 in both components of Euler's equation. From the r- 
component, equation (15), we obtain one relation between /o and g^ 2 : 

-pfg - Pcot9f' + /3 2 sin#cos# g'_ 2 + (2/3 2 cos 2 # - 2) g_ 2 = Ai + A 2 + A 3 ■ (27) 

Here, A±, A2, and .A3 are expressions involving fi and g-\\ 

, (/() 1 /l /l /r,o\ 

sin ft sm 9 

.1 cat - /! cot - /3 (/.J + /3 /f 0_! , (29) 

(30) 

From the ^-component, equation (16), we obtain a second relation between /q and (/_2 : 

+ P 5 L 2 - 2/3 2 sin# costf 9 _ 2 = B x + B 2 + B 3 , (31) 

where we have defined T> = 1 — /3 2 sin 2 0, and where the three terms on the righthand side are 
again combinations of f\ and g_y. 

_ 2 cotg fif[ 

°1 = /l ~ TH ~ 2^ ' ^) 

sm # sm 9 

£2 = Pfig-icate + pf[ g-i + ZPfig'-i , (33) 

B 3 = -2g_ x g'_ x . (34) 



Ax 


_ fi 


fifi<x»9 


sin 2 6 


sin 3 9 


A 2 


= Pfig- 


-i - 5 


A 3 


= 2 5 2 x - 


- 3ff_i . 



We have already found /i and g_i in the subsonic case of interest. After substituting these 
expressions, equations (25) and (26), into the righthand sides of equations (27) and (31), the coupled 
equations for /o and g~ 2 become: 



-Pfd - (3cot9f Q + p 2 sm9cos9 g'_ 2 + (2/3 2 cos 2 fl - 2) g_ 2 = 



+ 1 + y/V 



, (35) 



and 



-pfd + Vg'_ 2 -2/3 2 sin# cos9 <?_ 2 = (3 2 sm9cos9 



2 2 

P2 + ^3/2 



1 

v + 



1 



l + \/p) 2 



• (36) 



These last two relations constitute our second- order equations. For any value of /3, we may 
integrate them numerically from the upstream axis, 9 = ir, to the downstream axis at 9 = 0. 
Three initial conditions are required, of which we have already identified two: fo(n) = f'^ft) = 0- 
As a third initial condition, we use g- 2 (ir), whose value at this point is arbitrary. For each chosen 
value of g„ 2 (-7r), we may find fo(9) and g- 2 (9). We thus have a one-parameter family of outer flow 
solutions. 
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In the upper panel of Figure 3 we display, for the representative value (3 = 0.5, three solutions 
of g-2{0). We obtained each solution by assuming different values of <7_2(vr). Notice that g'_ 2 
vanishes on the upstream and downstream axes, implying again that the density profile is flat in 
both regions. Notice also that all curves attain the same value at 6 = tt/2. That is, g_2(7r/2) 
depends only on f3, and not on the prescribed initial condition 

The lower panel of Figure 3 shows the corresponding plots of fo(0). We see that /q(0) = 
in every case, ensuring regularity of u r on the downstream axis. This condition was not imposed 
a priori, but resulted automatically from integration of the governing equations. Specifically, the 
coefficient of /g in equation (35) includes cot 9, which diverges at = 0. Since all the other terms 
in this equation remain finite on the axis, /q(0) is forced to zero. 

Which of these solutions is the true outer flow for gas that is accreting steadily onto the 
gravitating mass? In principle, one could answer this question by continuing each solution inward, 
to see if the flow smoothly crosses the sonic surface, where u = 1. We shall not attempt such a 
calculation here. Instead, we will proceed by determining generically the mass accretion rate that 
is associated with each outer solution. Then, given the Bondi prescription for this rate, we will 
indeed be able to select the physical solution for each j3. 



4. Mass Accretion Rate 

4.1. Relation to Stream Function 

One could, in principle, equate coefficients of r~ 4 , r -5 , etc., and thereby obtain the coupled 
equations linking higher-order /- and ^-variables. We will now demonstrate, however, that the 
first- and second-order equations just presented are sufficient to establish the total accretion rate 
onto the mass. We will then relate, in Section 5 below, this infall rate to the desired friction force. 

Refer again to Figure 1 and imagine a sphere of radius r surrounding the mass. Reverting 
temporarily to dimensional variables, the mass accretion rate is 

M = -2vr / pu r r 2 sin 6 dO (37) 
Jo 

= 2, [V>(r,0) - i«r,7r)] , (39) 

where we have utilized equation (4) connecting u r and ip. Recall that ^(r, ir) is actually a constant, 
independent of r, and that we have set that constant to zero. We thus have 

M = 2vr^(r,0) . (40) 



To nondimensionalize this result, we first set the fiducial mass accretion rate to 2 7rpo c s r s- 
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Fig. 3. — Sample solutions to the second-order perturbation equations (35) and (36), for g_2{6) 
{upper panel) and fo(6) (lower panel). The initial values of g-2(it) are +1.12 (dotted), +0.12 (solid), 
and -0.88 (dashed). 
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After using the expansion of tp from equation (13), we obtain the nondimensional equation 

M = f 2 (0)r 2 + /i(0)r + / (0) + /-i^r" 1 + /_ 2 (0)r^ 2 ... . (41) 

One of our boundary conditions, ensuring regularity of uq on the downstream axis, is that fi(0) = 
for i = 1, —1, —2, etc. Since /2(C)) = 0, we find the simple relation 

M = /o(0) . (42) 

Both sides in this equation are functions of /?, although we have not indicated the dependence 
explicitly. In any case, the relation confirms our expectation that the mass accretion rate is inde- 
pendent of the sphere's radius r in steady-state motion. 2 We also now see that the higher-order 
variables /_2, etc. play no part in determining this rate. 



4.2. Relation to Density Perturbation 

Now that we have tied the mass accretion rate to /o(0), we can immediately rule out a subset of 
outer flow solutions as being unphysical. Figure 3 shows that, for g^ 2 {^) = 1-12, /o(0) is negative, 
corresponding to a net mass efflux. That such a situation is even possible emphasizes once more 
the need to extend the flow solution inward across the sonic surface. 

For this same choice of g- 2 (n), the dotted curve in the lower panel of Figure 3 shows that 
g~ 2 {0) < (7-2(71"). Indeed, we have just found one example of a general result: the difference 
<7_2 (0) — g~2( 7r ) agrees in sign with /o(0). We now show that the two quantities are in fact equal, 
apart from a multiplicative factor. 

Our proof starts with the fact that the lefthand side of the second-order equation (36) is a 
perfect derivative. Specifically, 

-Mo + Vg'_ 2 - 2 /3 2 sin cos 9 _2 = Jg(~Pfo + V 9-2) • (43) 

Turning to the righthand side of the same equation, we note first that sine? is an even function 
of 9 — vr/2, while cos# is an odd function. Since T> depends only on sin#, it has even symmetry. 
Inspection shows that the righthand side of equation (36) has odd symmetry. 

If we now integrate equation (36) from 9 = tt to 0, the righthand side vanishes because of the 
odd symmetry of the integrand. We find that 

C-yS/o + Vg. 2 ) d = ^ = (-Pf + Vg^ 2 ) e = . (44) 

Since /o(vr) = and T>(n) = T>(0) = 1, we have 

p_ 2 (7r) = -/9/o (0) + 3_ 2 (0) , (45) 



2 Note, however, that the original series expansion for i/j becomes inaccurate when r is not much greater than unity. 
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which we recast as 

/o(0) = ^ 2(0) = 9 ~ M • (46) 

Recalling equation (42) that identifies /o(0) as the mass accretion rate, we now see that this rate is 
proportional to the difference, upstream and downstream, of the second-order density perturbation. 
As f3 approaches zero, these two perturbations become equal. Indeed, the function g_2(#) is a 
constant (equal to 1/2) in the limit, consistent with a spherically symmetric flow. 



4.3. Modified Bondi Prescription 

To establish the physically relevant flow solutions, we need to specify the accretion rate as a 
function of velocity. Bondi (1952) fully solved the /3 = problem. That is, he determined the 
complete distribution of density and velocity surrounding a mass at rest within a background gas. 
Dimensionally, he found for the mass accretion rate 

c 

where A = e 3//2 /4 = 1.12 for the isothermal case of interest here. Recasting the rate into nondi- 
mensional form (recall Section 4.1), we have 



lim M = 2 A 

/3->0 



,3/2 

V > ( 48 ) 



as one constraint on the general form of M{(3). 

Prior to Bondi's work, Hoyle & Lyttleton (1939) studied accretion onto a mass traveling 
through a zero-temperature gas. Their dimensional result was 

4^0 G 2 M 2 , , 

M = — ^ . (49) 

Noting that the Mach number /3 is effectively infinite in this case, the equivalent, nondimensional 
finding is 



lim M = — . (50) 



2 

Bondi & Hoyle (1944) later showed that this relation provides an upper bound to the accretion rate 
in the zero-temperature case. Through more careful analysis of the wake, which here degenerates 
into an infinite-density spindle, they set the lower limit a factor of two smaller. 

The widely used interpolation formula of Bondi (1952) connects these limits, at least approxi- 
mately. Nondimensionally, the Bondi prescription is 

M(/3) = 1 — TN . (51) 

(1 + /3 2 ) 3/2 
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In the low-/? limit, M falls short of the isothermal result, but matches that for a 7 = 3/2 polytrope. 
The high-/? limit reproduces the lower bound established by Bondi & Hoyle (1944). 

Since we are focusing on the subsonic regime within an isothermal gas, we want our low-/? limit 
to agree with the exact result. Following Moeckel & Throop (2009), we adopt a modified form of 
the classic interpolation formula: 



1/2 

m w = " \: : : 2 / 2 > ( 52 ) 



2 ( A 2 + /? 2 ) 
(1 + /? 2 ) 



where we use the isothermal value of A previously given. For (3 <C 1, M approaches the result of 
Bondi (1952) given in equation (48). For /? > 1, we recover the upper limit of Bondi & Hoyle 
(1944). In the simulation of Moeckel &: Throop (2009) for an isothermal gas with /? = 10, this 
modified interpolation formula matches the calculated accretion rate to within 20 percent. Judging 
from their own polytropic simulations, both Hunt (1971) and Shima, Matsuda, Takeda, &; Sawada 
(1985) had earlier suggested that the original Bondi M(f3) be augmented by about a factor of two. 
In summary, equation (52) should be sufficiently accurate for our purposes. 

The combination of equations (42) and (52) gives us the proper value of /o(0) at each /3, and 
thus also establishes the physically relevant outer flow solutions. Figure 4 shows the physical /o(0) 
and <7_2(vr) as functions of /?. Note that the latter diverges as /? approaches unity. Thus, our 
perturbation series fails to describe the flow along the upstream axis in this limit. As we will show 
in the next section, however, the dynamical friction force remains finite for all j3. 

The three panels of Figure 5 display streamlines and isodensity contours for the indicated 
/?- values. These curves were constructed from equations (13) and (14) for i/j and p, respectively, 
using the three known terms in each series. The circle in each panel represents the sonic surface. 
As always, our results are only accurate well beyond this radius. 



5. Friction Force 

5.1. Integral Expression 

The dynamical friction force F is the total rate at which z-momentum is transferred from 
the background gas to the gravitating mass. Within our steady-state flow, the total momentum 
transfer rate into a surface surrounding the mass is independent of the size and shape of that 
surface, provided it lies outside the wake, where the physical interaction between the projectile and 
gas occurs. The net momentum flow calculated through such a surface integration all goes into the 
gravitating mass, causing its deceleration. 

Imagine the gravitating mass to be surrounded by a large sphere of radius r. In part, the z- 
momentum transfer arises from the advection of this quantity in the flowing gas across the spherical 
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Mach Number /? 

Fig. 4. — The upstream density perturbation g^i 7 ?) for the physical accretion flow, shown as a 
function of Mach number (3. This initial condition gives the correct M = /o(0), also shown in the 
figure. 
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Fig. 5. — Streamlines (solid) and density contours (dashed) for three different Mach numbers f3. 
The density contours correspond to p = 1.2, 1.4, and 1.6. The innermost streamlines enclose the full 
mass accretion rate M. Successive streamlines enclose 3, 5, and 7 times this rate. As in Figure 2, 
the inner circle represents the sonic surface. 
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surface. Since the inward flux of z- momentum is —pu r u z the kinetic portion of F is, dimensionally, 

rn 

F^[ n = —2tt / p u r u z r 2 sin 9 dO . (53) 



Another contribution to F is from the thermal pressure of the surrounding gas. This static portion 
of the force is 

/■7T 

Static = -2ir / p c 2 r 2 cos 9 sin 9 d9 . (54) 
Jo 

Adding these two pieces, we have, after nondimensionalization, 

rn rn 

F = — pu r u z r 2 sin 9 d9 — I pr 2 cos 9 sin 9 d9 , (55) 
Jo Jo 



where we have set the unit of force equal to 2 ir po c s rj. 

The integrand within the first, righthand term must be recast in terms of the stream function: 

pu u r 2 sm9 = (— V + — — — (56) 

r z pr 2 \d9 J pr 89 dr 

We may now evaluate F using the series expansions for ip and p. The full expression is a series of 
terms proportional to r 2 , r , r°, etc. 

All terms in F containing positive powers of r vanish upon integration. Those proportional to 
r 1 involve f\ and both of which are known explicitly. Terms associated with negative powers 
of r contain /- and (7-variables which we have not yet calculated (e.g., 5-3). However, as we 
consider ever larger radii r, where the series expansions themselves become increasingly accurate, 
these terms also go to zero. Only those independent of r survive. 

After restricting ourselves to r-independent terms, we find 

F = — [(1 - f3 2 ) sin0 cos 9 g_ 2 + P (l + cos 2 9) f' ]d9 . (57) 
Jo 

Here, we have omitted a number of terms in the integrand containing fi, g~\, and their derivatives. 
All of these terms are antisymmetric with respect to 9 — n/2 (i.e., they are odd functions), and 
therefore vanish upon integration. 



5.2. Relation to Mass Accretion Rate 

By dimensional considerations, the friction force should be F = CMV, where C is dimen- 
sionless. In the hypersonic limit, this multiplicative factor contains a Coulomb logarithm (e.g. 
Ruderman & Spiegel 1971). We now demonstrate the surprising fact that, in the subsonic case of 
interest here, the factor is exactly unity. In fully nondimensional language, we shall prove that 



F = M/3 



(58) 
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We begin by splitting the integral on the righthand side of equation (57) into two parts: 

/*7T f'TT 

F = — I PfodO- / [(1 -/3 2 ) sin0 cos 9 g_ 2 + /3 cos 2 9 fd]dd (59a) 
Jo Jo 

= /3/ (0)-X (59b) 
= PM-I. (59c) 

In these equations, we have used the fact that fo(ir) = and /o(0) = M (eq. 42). We have further 
defined 

r-7T 

T = / d9 [(1-/3 2 ) sin^ cos 9 g- 2 + P cos 2 6 f ] . (60) 
Jo 

We next show that I vanishes. 

First recall that our flow is irrotational. Specifically, the (^-component of the vorticity vanishes, 
so that 

a0 or 

Expressing both velocity components in terms of the stream function through equations (4) and 
(5), we have 

89 d9 + pCOtt> 89 9 89i + r dr dr P r dr 2 " u • 

We substitute the series expansions for ip and p into this last equation and set the coefficients 
of all powers of r to zero. Following this procedure for r 2 and r , and using the known expressions 
for f 2 , fi, and <?_i, yields identities. However, setting the r-independent terms to zero leads to a 
nontrivial result: 

1 \f 

Pfg - cote f' Q - p 2 sin 9 cos 9 g'_ 2 + 2 p 2 sin 2 6 g- 2 = ~ p • (63) 

We add this last equation to the second-order equation (35), obtaining 

2 (P 2 - 1) g_ 2 - 2/3cote/^ = | - JL + . (64) 



(65) 



Multiplying through by —(1/2) sin e cos gives 

/ 1 9 

(1 - p 2 ) sine cos e 5 _ 2 + Pcos 2 ef Q = -sinecose^- - — + 



1 + v 7 ^ 



Integrating over e, we recognize the lefthand side of the resulting equation as I. The righthand 
side vanishes, since the integrand is an odd function. We see therefore that equation (58) holds. 

If we now employ the modified Bondi prescription, equation (52) for M, we have an explicit 
expression for the force: 

2(3 (A 2 + p 2 ) 1/2 , 

F = -4 . 66 

(1 + P 2 f 
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Figure 6 displays the function F(/3). Also shown, as the dashed curve, is the result from Ostriker 
(1999) in which the force diverges as f3 approaches unity. In the limit of low /3, both forces rise 
linearly from zero at /3 = 0, but our initial slope is larger by a factor of 3A = 3.36. Indeed, over 
most /3-values, our force exceeds that derived by Ostriker (1999), presumably because we have 
included both the gravitational tug from the wake and the direct accretion of momentum from the 
flow. Our force does not rise monotonically but instead peaks around /3 = 0.68 and then begins to 
decline; we expect this decline to continue into the supersonic regime. We should bear in mind that, 
while equation (58) is exact, equation (66) for F is only as accurate as the underlying interpolation 
formula. 



6. Velocity and Mass Evolution 

Our simple result for the dynamical friction force means that the deceleration of the gravitating 
mass is also simply described, as long as there are no other forces at play. As we have stressed, the 
force is the rate at which gas transfers linear momentum to the object. But the object's momentum 
is MT^, where we now revert to dimensional variables. In the reference frame where the background 
gas is stationary, we have 

« = -MV, (67) 

which implies that 

1 dV 2 dM 

V~dt ~ ~M ~dT ' ^ ' 

If Vq and Mq are the object's initial speed and mass, respectively, then 



V ( M 



- 



-2 



(69) 



Vq \M , 

To track the speed as a function of time, we rewrite equation (52) for the mass accretion rate 

as 

dM 4irp c s r 2 s (A 2 + /3 2 ) 1/2 / M 



dt (1 + /3 2 ) 2 \M J ' (70) 

Here, f3 = V/c s as before, while r s is now defined in terms of the initial mass: r s = 2G Mq/c 2 s . 
The fully nondimensional evolutionary equation for the speed is then 

1\ df3 _ 4(A 2 + /3 2 ) 1/2 ( p\- 1 ' 2 



Pj dr (1 + £2)2 V/V 

In this last equation, we have introduced the initial, nondimensional speed /3o, as well as a nondi- 
mensional time, t = t/tQ, where 

*° = o %2 i\/r ( 72 ) 

2 7r po G z Mq 

= H 5 *-*- < 73 > 

2 vr po c s rj 
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Fig. 6. — The dimensionless friction force F as a function of Mach number /3. The dashed curve 
shows the force derived by Ostriker (1999), which diverges as j3 approaches unity. 
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The denominator in equation (73) is the fiducial mass accretion rate defined in Section 4.1. Thus, 
to is of order the accretion time onto the initial mass. 

The upper panel of Figure 7 plots /3(r) for /3q = 0.2, 0.5 and 0.8, obtained by numerical 
integration of equation (71). Also shown, in the lower panel, is the growth of the nondimensional 
quantity M, the mass of the gravitating body relative to its initial value. As expected, the body 
slows down appreciably within an accretion time. 



7. Summary and Discussion 

This study has pivoted on the close relationship between the dynamical friction force, i.e., the 
transfer of linear momentum from gas to a gravitating object, and the transfer of mass to that same 
object. This relationship is embodied in our central result, equation (58). From this equation, in 
turn, we derived an analytic expression for the force itself, equation (66). 

We are now in a position to address a basic question raised in Section 2.1. How are we justified 
in assuming steady-state flow, when the gravitating body is continually decelerating? The answer 
is that quasi-steady flow is established within a radius r cr ; t over which the sound crossing time 
(^crit/cs) equals the time for the object's momentum to decrease appreciably (MV/F). Recalling 
that F is normalized to 2tt poc^r^ and using equation (58), we have, nondimensionally, 

r C rit = ol , (74) 

where a = Mq/(2 tt po r^). The latter quantity was implicitly assumed to be large from the start, 
when we neglected the self-gravity of the gas. The nondimensional mass accretion rate M = /o(0) 
hovers near unity for the entire evolution (recall Fig. 4), while M itself starts at unity and climbs. 
Hence, the critical radius is much larger than r s , and our analysis is self-consistent. 

We note that dynamical friction still operates in circumstances where mass accretion is frus- 
trated. For example, a wind-emitting star moving through a gas cloud experiences mass loss rather 
than mass gain. Cloud gas impacting the wind upstream is arrested or refracted in a bowshock, as 
analytically calculated by Wilkin (1996). Downstream, the wind forms a supersonic jet. As long 
as the upstream standoff radius of the shock lies within r s and the downstream jet is relatively 
narrow, the far-field perturbations are close to what we have obtained, and equation (66) for F still 
applies. 

When the object is actually able to accept gas freely, dynamical friction arises in two physically 
distinct ways. First, there is the gravitational tug from the wake. Second, momentum is transferred 
directly to the object by gas falling onto it. Our finding that these two forces sum to M V is at least 
roughly consistent with simulations. In a numerical study directed primarily at the mass accretion 
issue, Ruffert (1996) explicitly determined both force contributions on accretors of various size 
in a 7 = 1.01 gas. For R/r acc = 0.1 and (3 = 0.6, the simulation ended before the flow reached 
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Fig. 7. — Evolution of a particle's speed and mass as a function of nondimensional time r. The 
different curves represent initial speeds j3o = 0.8, 0.5, and 0.2. A particle both triples its mass and 
slows to ~ 0.1 times its initial speed in a fraction of its mass accretion time. 
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steady-state (see his Fig. 2). After initial transients died out, the gravitational drag was steady 
until t ~ 13 iBH> where the Bondi-Hoyle time £bh 

is ?" acc /c s . Thereafter, this force component 
declined for the rest of the integration. At the end of the simulation (t = 32 £bh)> the sum of 
the gravitational drag and momentum accretion forces was 1.2 MV. For R/r acc = 0.02 and the 
same Mach number, the two forces quickly leveled off, with a sum equal to 1.4 M V . However, this 
simulation ran only until t = 10 iBH> so it is not clear whether the gravitational drag would have 
later declined, as in the first case. 

Following historical precedent, we have restricted our investigation to an isothermal gas. For 
an isentropic gas with 7 > 1, it seems likely that the friction force will still be given by MV, as 
long as the accretor is moving subsonically. Verifying this equality analytically would require a 
perturbation study analogous to the present one. We leave such a project for future investigators. 

Again the current body of numerical studies is in broad accord with our expectation. Ruffert 
(1994) determined the total friction force on an accretor moving through a 7 = 5/3 gas. For 
R/r^cc = 0.1 and (3 = 0.6, the friction force was 1.1 MV at t = 70 £bh- For R/r acc = 0.02 and 
the same Mach number, the flow had not achieved steady state by t = 19 £bh- The total force 
was 1.8 MV at this time, but was falling rapidly. A future project of interest would be to redo 
these simulations over a range of (3- and 7- values, running the simulations long enough until a true 
steady state is reached. 

For the more general isentropic case, M can no longer be approximated by equation (52). 
Instead the value of M at a given V decreases with higher 7- values, as shown analytically by Bondi 
(1952) for V = 0, and as seen in the simulations of Ruffert (1994, 1995, 1996) for accretors moving 
relative to the background gas. Isentropic flows are less compressible than isothermal ones, so the 
wake will be less dense. As a result, the friction force will also be lower, presumably by the same 
amount as the accretion rate M. 

In the present investigation, we have been unable to tease apart analytically the two force 
contributions. To do so would require study of the flow closer to the gravitating mass, specifically 
across the sonic surface. In principle, a perturbation series in this region could be linked to the 
outer one developed here. Besides elucidating the momentum transfer through infall, such a study 
could also establish M analytically as a function of velocity, thus putting accretion theory as a 
whole on a firmer foundation. 
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